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We describe a variational approach to solving Anderson impurity models by means of exact di¬ 
agonalization. Optimized parameters of a discretized auxiliary model are obtained on the basis 
of the Peierls-Feynman-Bogoliubov principle. Thereby, the variational approach resolves ambigui¬ 
ties related with the bath discretization, which is generally necessary to make Anderson impurity 
models tractable by exact diagonalization. The choice of variational degrees of freedom made here 
allows systematic improvements of total energies over mean field decouplings like Hartree-Fock. 
Furthermore, our approach allows us to embed arbitrary bath discretization schemes in total energy 
calculations and to systematically optimize and improve on traditional routes to the discretization 
problem such as fitting of hybridization functions on Matsubara frequencies. Benchmarks in terms 
of a single orbital Anderson model demonstrate that the variational exact diagonalization method 
accurately reproduces free energies as well as several single- and two-particle observables obtained 
from an exact solution. Finally, we demonstrate the applicability of the variational exact diagonal¬ 
ization approach to realistic five orbital problems with the example system of Co impurities in bulk 
Cu and compare to continuous-time Monte Carlo calculations. The accuracy of established bath 
discretization schemes is assessed in the framework of the variational approach introduced here. 

PACS numbers: 72.80.Rj; 73.20.Hb; 73.61.Wp 


I. INTRODUCTION 

The Anderson impurity modeli (AIM) is a general 
model for the description of interacting impurities in 
metallic host systems. Originally, it was developed to de¬ 
scribe single atoms with open d- or f-shells embedded in 
bulk materials and to understand the formation of their 
magnetic moments^. Furthermore the model includes the 
widely discussed Kondo physics^. Multi orbital variants 
of the AIM gained considerable attention in the context 
of rare-earth impurity systems^ as well as more recently 
magnetic adatoms or molecules on surfaces^—. Finally, 
dynamical mean held theory^ (DMFT) links correlated 
bulk systems as well as nanostructures to Anderson im¬ 
purity models. 

To address the electronic structure of realistic corre¬ 
lated electron materials one often resorts to LDA+-1- 
approaches^, where quantum lattice or impurity mod¬ 
els are derived from first principles calculations. The re¬ 
sulting models are typically multi-orbital models includ¬ 
ing complex hybridization between the impurity and a 
continuous bath of states from the surrounding material, 
which brings along two challenges: First, the numerical 
solution of the impurity models and second the interpre¬ 
tation of the physics contained in these generally complex 
models in more simple terms. Experiments are for in¬ 
stance often interpreted in terms of atomic spins, crystal 
Held, ligand fields or cluster approaches^, which typi¬ 
cally involve a small discrete set of bath states or no bath 
states at all. The link of the complex, ab initio derived 
models and simpler phenomenological models is a priori 
unclear and relates to the so-called bath discretization 
problem of exact diagonalization solvers of the AIM. 


The solution of the Anderson impurity model for gen¬ 
eral parameters has to be done numerically by means of 
e.g. quantum Monte Carlo^ (QMC), numerical renor¬ 
malization groupii (NRG), or exact diagonalization (ED) 
methods^ While NRG and QMC are in principle nu¬ 
merically exact methods, they become computationally 
very demanding, when dealing with many orbitals, hy¬ 
bridization functions with low symmetry, spin orbit cou¬ 
pling and general fermionic four operator Coulomb ver¬ 
tices. ED methods deal with low symmetries and general 
Coulomb vertices at no additional computational cost 
but suffer from the so-called bath discretization prob¬ 
lem: Due to the exponential growth of the many particle 
Fock space with the system size, it can handle only a 
few bath levels per orbital. A mapping of the continu¬ 
ous bath to a discrete version has to be found. Several 
approaches to this task have been introduced. One is to 
fit the hybridization function of the continuous bath on 
Matsubara frequencies^, another is to represent the hy¬ 
bridization function by a continued fraction and to link 
its coefficients to the parameters of the bathi£. These 
schemes are systematic in the sense that they converge 
to the full model when including more and more bath 
sites. However, in the multi-orbital case, the number of 
bath sites is limited (typically on the order of three or 
less for a five orbital impurity problem), so the quality of 
the mapping can hardly be checked by an analysis of the 
convergence. 

Basically two different strategies have been laid out to 
circumvent this problem. First, the many body Hilbert 
space can be truncated in the sense of configuration in¬ 
teraction (Cl) expansions, which have a long tradition 
in the context of quantum impurity problems^ and are 
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subject of recent developments— d Cl expansions are 
variational, i.e. they deliver upper bounds for total ener¬ 
gies, but they do not provide simplified auxiliary Hamil¬ 
tonians. On the other hand, there are several approaches 
towards optimized cluster approximations to Anderson 
impurity problems. In this context, self-energy functional 
theory2£ is based on an extremal principle but it is not 
variational regarding total energies and does not allow 
for variations of interaction parameters or the interacting 
orbitals. More general optimizations are possible in the 
framework of the so-called self-energy embedding theory 
(SEET)2I, which is however not variational. 

In this paper, we combine ideas of variational ap¬ 
proaches and optimized cluster approximations to the 
AIM. We introduce a strictly variational method of ap¬ 
proximating an AIM with continuous bath by an AIM 
with finite strongly reduced number of bath sites, which 
we call variational ED method. It guaranties an optimal 
approximation to the AIM for a given number of bath 
sites in the sense of thermodynamic ground state proper¬ 
ties. The method is based on the well-known Peierls- 
Feynman-Bogoliubov variational principle^—, which 
finds optimal effective models on the basis of an optimal 
density matrix by minimizing a free energy functional. 

We will introduce the AIM and the variational prin¬ 
ciple in Sec. eh where we also explain the details of 
calculating the Peierls-Feynmann-Bogoluibov free energy 
functional and how to minimize it efficiently. By treating 
a single orbital model with the variational ED method 
in Sec. m we analyze its performance in comparison 
to an exact treatment, established bath discretization 
methods^ as well as Hartree-Fock theory. In Sec. m 
we demonstrate the applicability of the method to realis¬ 
tic five orbital system with the example of Co impurities 
in bulk Cu and compare to QMC simulations. We show 
that the variational ED method leads to systematically 
lower, i.e. more accurate, free energy estimates than un¬ 
restricted Hartree-Fock and traditional bath discretiza¬ 
tion schemes also in the multiorbtial case. Conclusions 
and outlook are given in Sec. IVII 


II. MODEL AND METHOD 


After introducing the Anderson impurity model we will 
recapitulate the Peierls-Feynman-Bogoliubov variational 
principle and show how to apply it to discretize Anderson 
impurity models in an optimal manner. 


A. The Anderson impurity model 


The bath is described by 

-Hbath = 'y ' ^akl^aka^ (2) 

Oik,<7 

where e ak is the energy of the bath state with 
band/orbital index a and some additional quantum num¬ 
ber k. n c ak(7 = c\ xk(J c ol ka is the corresponding particle 
number operator. The hybridization part 

Hhyb = y ' bafc (3) 

Oik,C 7 

couples the bath sites of one band to an orbital of the 
impurity with a coupling strength I 7 ak . The bath elec¬ 
trons with spin a are created and annihilated by c^ fco . 
and c a k a , respectively, while (fi aG ( d aa ) denote the cre¬ 
ation (annihilation) operators of the impurity electrons. 
The impurity site is described by 

Hjmp — y ' £a n aa T ^ ' U a pry^d^ aa dg a / d^c' dga , 

ct,a oi,fi , 7 ,8,<7<j' 

(4) 


which contains the on-site Coulomb interaction U a pjg 
and the on-site energies E d a . By integrating out all bath 
degrees of freedom we arrive at the hybridization function 


A a (ui) 
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v: k v ak 
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(5) 


which describes the energy dependent coupling of the im¬ 
purity to the bath. 


B. Peierls-Feynman-Bogoliubov variational 
principle 

Given a Hamiltonian H , which is “difficult” to solve, 
we search for an optimal approximation to H within a 
set of simpler effective Hamiltonians H. The Peierls- 
Feynman-Bogoliubov variational principle^— provides 
us with a prescription on how to fix the parameters of H 
in a thermodynamically optimal way, i.e. such that the 
canonical density matrix resulting form H approximates 
the density matrix corresponding to H as close as possi¬ 
ble. More strictly speaking: the canonical density opera¬ 
tor pjj = 1/Zjj exp(— /3H) of the auxiliary system, where 
Zjj = Trexp(— /3H) is the partition function, approxi¬ 
mates the exact density operator p derived from H as 
close as possible, when the Peierls-Bogoliubov-Feynman 
functional 


The Hamiltonian of the initial AIM (termed “original 
model” hereafter) reads 


H — Hbath T Hhyb T 


( 1 ) 


*\Ph\=*h + (H-H)h, (6) 

becomes minimal. Here = —-jin Zpj is the free en¬ 
ergy of the effective system. (H — H )^ = Tr p^(H — H) 
denotes a thermodynamic expectation value with respect 



3 



A 


) ^ effective model 



Figure 1. (Color online) Illustration of the original and effec¬ 
tive model for the case of one orbital and six bath sites. Blue 
represents bath character and red impurity character: In the 
effective model bath and impurity states can be mixed. e)) nc 
are eigenvalues of h£ k ,. 


to the effective system. In the case of Pfi = p the func¬ 
tional becomes minimal and coincides with the free 
energy of the original system. In our case H repre¬ 
sents the full AIM, Eqs. (|T ]) - (p[j) . and H is the model with 
discretized bath, which is now introduced. 


C. Effective Hamiltonian 


The structure of the effective Hamiltonian for the case 
of a single impurity orbital is depicted in the right panel 
of Fig. [I] In contrast to the original model (left panel 
of Fig. [T]). the effective model consists of two decoupled 
parts: First, the effective impurity coupled to one bath 
site only and second the remaining bath sites. I.e. we 
partition the full Hilbert space T~L into a correlated sub¬ 
space C (first part) and an uncorrelated rest TZ (second 
part). In this work, we consider for concreteness a cluster 
consisting of a multi-orbital impurity and one bath site 
per impurity orbital for the correlated space but other 
choices are similarly possible. The single particle states 
of the effective model are related to those of the origi¬ 
nal model by a unitary transformation, which allows for 
mixing of original “bath” and “impurity” character in 
the effective model. 

The optimal matrix elements of the effective model, as 
well as the optimal unitary transformation are found by 
minimizing the functional j6]). The states spanning C are 
defined by 


K) = U d dl\d a ) ( 7 ) 

k 

|cal) = U C d a J \d a ) + Y |Cafc>, (8) 

k 

where the coefficients u are chosen such that \d a ) and 
|c Q i) form an orthonormal basis of C. An orthonormal 


basis spanning TZ is defined by 

| c ak ) = u c d “ k \d a ) + Y u cZ k ,,\ c <xk'),k > 1. (9) 

k' 

As a whole, the coefficients u form a unitary matrix. In 
practice, we obtain the elements of this matrix from the 
QR decomposition of a matrix, in which the first two rows 
are defined by the coefficients of \d a ) and |c a i) and all 
other elements are zero. This leads to a new orthonormal 
basis for the full space TL 1 which provides the partitioning 
according to TL = C © 1Z. The ansatz for the effective 
Hamiltonian in this new basis explicitly reads 

H = H C + H n , (10) 

with 

H — ^ ^ hh ^d a ^ ^ ^ctiCgiCcti 

a _ _ . a . (11) 
~h ^ v ^a.dgdg H~ ^ ^ U 0 t.p^sd\ yL(J dp a ,d 7<t' dga 
cx. 

and 

H n = Y h «kkAk^k'- ( 12 ) 

cx,(k,k')> 1 

It is stressed, that the new states are linear combi¬ 
nations of the original impurity and bath states, lead¬ 
ing to mixed basis states. The new impurity states can 
have some amount of bath character and vice versa. The 
Hamiltonian in Eq. (EH) states a many-body problem 
which can be solved by exact diagonalization, as long as 
its Hilbert space is sufficiently small. In contrast, the 
Hamiltonian m states a one-particle problem and can 
be solved by diagonalizing the matrix h^ kk ,. In sum¬ 
mary, the Hamiltonian H = H c + H 71 defines an effective 
Hamiltonian, which can be solved exactly and thus the 
functional ([6| can be calculated. 

This ansatz implies several approximations. First, all 
couplings between C and 7 Z are neglected. Second inter¬ 
action terms are restricted to new effective impurity or¬ 
bitals d a within C, which is motivated by the fact that the 
original model includes only on-site interactions too. The 
latter approximation can be relaxed to include arbitrary 
interactions within C, but we keep it here for simplicity. 

Finally, we note that the amount of variational degrees 
of freedom in the variational ED approach is such that 
it includes Hartree-Fock as the limiting case Uafi-yS —>• 0. 
Thus, we expect that variational ED will generally give 
more accurate energy estimates than Hartree-Fock. 


D. Implementation 

In order to perform the minimization in practice, the 
number of free parameters has to be kept sufficiently low. 
First, for the rest of this work it is assumed that the 












4 


Coulomb tensor U ai 3-yS is not varied. We choose it to be 
the same as in the original model. Test calculations have 
shown, that the variation of the Coulomb tensor is not 
crucial, as this can mostly be absorbed into the variation 
of the impurity level. The single particle matrix elements 
of H c are assumed to be free parameters. In principle, 
the parameters of the uncorrelated Hamiltonian are free 
parameters, too. However, to further reduce the number 
of free parameters, we define h^ kk , by a projection of 
a Hartree-Fock solution of the original Hamiltonian onto 
the states \c a k)- The Hartree-Fock solution of the original 
Hamiltonian (0 can be written as 

#HF = J2 £ n F( icn, (13) 

n 

where the eigenstates |n) and energies e^ F are found by 
applying the Hartree-Fock decoupling 

— {dacrd"y<T')d'p lT ids a — (dg cr /d,5 (7 }d£, 0 .cLy 0 .' 

(14) 

to 0 and solving the resulting non-interacting prob¬ 
lem self-consistently. The single particle matrix elements 
within the uncorrelated space 1Z explicitly read 

h akk'=J2^n F {c a k\n)(n\c ak '). (15) 

n 

In order to not break any spin rotation symmetries, re¬ 
stricted Hartree-Fock is used. 

The functional 4>[p^] now depends on the unitary 
transformation and on the matrix elements of H c . The 
minimum of the functional is searched by iterative meth¬ 
ods. Thus, the functional l>[p^] has to be calculated for 
various points of the variational space with the computa¬ 
tionally most expensive part being here the diagonaliza- 
tions of H c . Therefore, we first search for fixed parame¬ 
ters in H c a corresponding optimal unitary transforma¬ 
tion matrix defining the optimal partitioning T~L = C © 1Z 
using an SLSQP algorithm^. The search of the mini¬ 
mum w.r.t. the parameters of H c is then done by the 
Nelder-Mead algorithm^. The number of independent 
parameters can be further reduced when the original sys¬ 
tem shows symmetries like orbital degeneracies which are 
assumed not to be broken in the effective model. 


III. BENCHMARK FOR A SINGLE ORBITAL 
AIM 

In this section the variational ED method is tested for 
its performance in reproducing the density operator as 
well as observables such as the occupation number, dou¬ 
ble occupancy and crystal orbital overlap populations of 
a simple original model. The original model we consider 
here is a single orbital model with only 6 bath sites, which 


itself can be solved by exact diagonalization. The de¬ 
tailed setup of the model is as follows: The impurity level 
is £d = —2.0 eV, the interaction strength is U = 4.0 eV. 
The 6 bath levels are equally aligned around a mean bath 
energy £& in an interval of 2 eV (i.e. the bandwidth of the 
bath). The coupling is 14 = 0.9 eV. The mean bath en¬ 
ergy £b is swept from —6.0 eV to 6.0 eV. All energies are 
measured w.r.t. to the Fermi energy £p = 0. The system 
is solved for T = 0. The model is first solved exactly, sec¬ 
ond by the variational ED method, third by unrestricted 
Hartree-Fock and finally by ED using reduced bath sites 
obtained by fitting of the hybridization functions on the 
imaginary Matsubara frequency axisi^. The latter type 
of approaches require generally the introduction of a so- 
called weight function W n for the fitting procedure, as 
explained in the appendix ©■ 

The central object for the assessment of the quality 
of the methods is the difference between and the 

exact free energy $#, as shown in Fig. 0). For bath 
sites energetically far away from the Fermi level and from 
single particle excitation energies of the impurity (|e&| > 
4 eV), all methods lead essentially to the correct free 
energy. Deviations occur, however, for bath levels closer 
to the Fermi energy. The Hartree-Fock free energy differs 
from the exact thermodynamical potential on the order 
of 100 meV basically in the whole range of |et,| < 4 eV. 
The fitting of the hybridization function function on the 
imaginary axis leads to rather accurate free energies as 
long as all bath sites are above or below the Fermi energy 
(|£b| > 1 eV), while for |£&| < 1 eV deviations from the 
exact thermodynamical potential on the order of 10 to 
20 meV occur. The choice of an optimal weight function 
(see appendix E| depends on details of the bath: For 
the case of bath sites on both sides of the Fermi level 
(|£b| < 1 eV) W n = l/w„ leads to the lowest free energies. 
Otherwise, the constant weight function W n = 1 shows 
smallest deviations of the free energy functional from the 
exact solution. The variational ED method is generally 
very close to the exact solution. Only for the special 
case of a strictly symmetric distribution of the bath sites 
around the Fermi energy (£& = 0 eV) a deviation on the 
order of meV occurs. 

Fig. Hb)-e) shows a comparison of several observables 
(chemical bond strength b), bath occupation c), impurity 
occupation d) and double occupation e)) calculated with 
the different methods. For the outer most regions (|£&| > 
4 eV) all methods describe the observables accurately. 
The fitting of hybridization functions on the Matsubara 
axis leads to deviations depending on the weight function, 
especially for the double occupation and chemical bond 
strength if the bath is centered around the Fermi energy 
(£f, ss 0 eV). Hartree-Fock systematically overestimates 
the double occupancy for |£t,| <4 eV. The variational 
ED method shows nearly no deviations from the exact 
solution at all. 

It is instructive to examine the unitary transformation 
linking the basis of the original and effective model. Fig. 
[2f) shows the coefficients of the linear combination of the 
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Figure 2. (Color online) Benchmark of different ED ap¬ 
proaches and spin-polarized Hartree-Fock theory against an 
exact solution for single orbital Anderson impurity models 
with a mean bath energy £(,. (a) Difference between the free 
energy functional obtained by different approximate meth¬ 
ods according to Eq. and the free energy of the original 
model. The inset shows a close up view for the vicinity of 
the Fermi energy on a logarithmic scale, (b)-(e): Compari¬ 
son of local and non-local observables obtained from an exact 
solution (“orig”, bold cyan) and calculated by the four dif¬ 
ferent approximate methods, i.e. the variational ED method 
(“var”, solid black), Hartree-Fock (“HF”, dashed red) and fits 
of hybridization functions on the imaginary axis with different 
weight functions (W n = 1: “fitO”, dotted blue; W„ = l/ui n : 
“fitl”, dashed blue.) Panel b) shows the chemical bond 
strength, c) shows the total bath occupation, d) the impu¬ 
rity occupation and e) the double occupation. Panel e) shows 
the coefficients of the unitary transformation linking the orig¬ 
inal model to the optimized effective model with one bath-site 
per spin-orbital in C for the example of Eb = —0.3 eV. 


states spanning the correlated space | d) and |£i) (see Eqs. 
(0 and ®) for an original model with the bath centered 
around the energy et, = —0.3 eV. The effective impu¬ 
rity has mainly | d) character with small bath admixture 
and can approximately be interpreted as the old impurity 
state. The coupled effective bath state is nearly a pure 
linear combination of old bath states, where states closer 
to the Fermi energy contribute stronger than those fur¬ 
ther away. This behavior is very reminiscent of effective 
bath wave functions obtained in variational approaches 
like the Varma-Yafet^ or the Gunnarsson-Schonhammer 
expansion! 4 

For the treatment of original models with far more bath 
sites, it is important to note, that the coefficients defining 
the unitary transformation from the original bath states 


to the effective impurity and bath orbitals, i.e. u^} and 
vary smoothly as function of the bath energies on 
either side of the Fermi energy. 


IV. CO IMPURITIES IN CU: APPLICATION 
TO A REALISTIC FIVE ORBITAL SYSTEM 

A. The AIM derived from LDA 


In this section the variational ED method is applied to 
a realistic model of Co impurities in bulk Cu, which has 
been obtained from super-cell DFT calculations and has 
been analyzed using a QMC impurity solver in Ref. @. 

The cubic symmetry of the Cu crystal leads to a split¬ 
ting of the Co 3d-orbitals into blocks of t 2 g and e g sym¬ 
metry. From the DFT hybridization function, which is 
a continuous function, we obtain our initial model as¬ 
suming some large number of bath sites, here 100 per 
orbital. (This number does not present a limiting fac¬ 
tor and could be chosen arbitrarily larger). The bath 
sites are assumed to be equidistantly distributed between 
— 10 eV and 10 eV, and the hybridization terms V a k are 
then found by fitting the imaginary part of a discretized 
hybridization function 


^disc(^) — ^ ^ 
k 


VakVak 
u - £k + iS' 


(16) 


with some broadening S = 0.1 eV to the ab initio hy¬ 
bridization function A(w) on the real axis. The VA are 
plotted in Figurc[3ji). The crystal field obtained from the 
DFT calculation is £g g — e/ 2g = 0.136 eV. As in Ref. 3, 
we consider a rotationally invariant Coulomb interaction 
defined by 


21 

Uaf_ 3-y<5 — y ^ Ufc (Om/^m, ^)m$rn)F ? (17) 

fc=0 


where afc(a m /3 m ,7 m i5 m ) are the Gaunt coefficient s 27 ! 28 
and where F° = U, F 2 = 14/(1 + 0.625)J and F 4 = 
0.625 F 2 are Slater parameters with the average Coulomb 
interaction U = 4.0 eV and Hunds exchange interaction 
J = 0.9 eV. Due to the so-called double counting prob¬ 
lem inherent to LDA+-(- approaches, the filling of the im¬ 
purity d-levels is not exactly known. Here, we consider 
the double counting potential [i = 27 eV as in Ref. 3- 
All data is obtained at a inverse temperature of /3 = 40, 
like in the case of the QMC simulations. Finally, we as¬ 
sume that the cubic symmetry of the system prevails, 
which means that only two independent sets of matrix 
elements (for the t 2 g and e g states) have to be varied 
during the minimization of 































































6 


B. Implementation of the variational ED method 
for the 5 orbital AIM 

We compare two different sets of variational degrees 
of freedom for the optimization of the one particle basis, 
which we refer to as “bath” and “all”. In the “bath” case, 
only bath sites are optimized, i.e. we fix the expansion 
coefficients it^ 1 = 0, = 0 and = 1. This leads 

to considerably less variational parameters and a much 
smaller amount of expectation values to be calculated in 
each step of the iteration. In the second approach, “all”, 
which is computationally more demanding because the 
full two-particle density matrix of the effective system 
has to be calculated, we optimize the full one particle 
basis of the bath and that of the impurity. 

Because a full optimization of the parameters of the 
effective model is computationally challenging, it is cru¬ 
cial to start the optimization from a good initial guess. 
We obtain such initial guesses for the parameters of the 
bath by fitting of hybridization functions on Matsubara 
frequencies as introduced in the appendix lAl We choose 
= £ a as the initial guess for the parameters of the im¬ 
purity. The resulting first guesses using different weight 
functions are summarized in the Table |U While all weight 
functions lead to setups with the t 2 g bath sites above the 
Fermi energy and the e g bath sites below, the details of 
their energetic positions and the hybridization strengths 
depend strongly on the form of W n . Adding more weight 
on features on small Matsubara frequencies shifts the ef¬ 
fective bath parameters to smaller values. The quality 
of these starting guesses in the context of the variational 
principle is discussed in the next section. 

Table I. Parameters of the effective model (see Eq. CD) 
obtained by the fit of hybridization functions on imagi¬ 
nary frequencies using different weight functions (W„ = 
1,1 /uJn, 1/4, see appendix [All and the iterative optimization 
(“var”). 


weight function 

1 

1 /UJn 

1 /Wn 

var 

et 2g i(eV) 

3.203 

0.775 

0.068 

2.658 

Ut 2g (eV) 

1.563 

0.606 

0.223 

1.717 

£e g i(eV) 

-2.314 

-0.019 

-0.015 

-1.995 

K g (eV) 

1.049 

0.170 

0.156 

1.418 

4 g (eV) 

-27.30 

-27.30 

-27.30 

-26.57 

4(eV) 

-27.44 

-27.44 

-27.44 

-27.87 


The large number of bath sites (100 per orbital) in the 
original model leads to 202 variational parameters defin¬ 
ing the unitary transformation in the “all” case or 100 
parameters in the “bath” case for each orbital. The ob¬ 
servation, that u d c 2 k , are smooth functions of energy 
(c.f. Fig. [2t) below and above Ep, leads to the possibil¬ 
ity of expanding them in a set of smooth functions and 
thereby reducing the number of variational parameters 
considerably. Here, we chose five Chebyshev polynomials 
T n (k) per orbital for bath sites above and five for those 


Table II. The free energy functional <f>, the total impurity oc¬ 
cupancy rid. and the local spin S as obtained from simulations 
of the AIM for Co impurities in Cu. The values of <f> are 
shown as differences to the results from unrestricted Hartree- 
Fock (UHF): AT = 4> — $uhf- Total impurity occupation and 
spin calculated with the variational ED method are compared 
to QMC solutions of the AIM from Ref. @. Different flavors of 
the variational ED method are considered: first “bath” and 
second “all” with the model parameters obtained from the 
fits of the hybridization function on the imaginary frequencies 
using different weight functions ( W„ = 1 “fit0”,Wn = l/w„ 
“fitl” and W„ = 1/hJn “fit2”) and finally full optimization 
of transformation and model parameters labeled “var”. Re¬ 
stricted Hartree-Fock (RHF) and unrestricted HF (UHF) re¬ 
sults are also shown. 



Af-(eV) 

(rid) 

S 

QMC 


7.78 ± 0.05 

0.92 ± 0.02 

UHF 

0 

7.78 

1.78 

RHF 

0.52 

8.20 

1.06 

fit0,bath 

0.05 

7.75 

1.06 

fitl,bath 

1.15 

7.84 

1.04 

fit2,bath 

2.77 

7.92 

1.03 

fit0,all 

-0.22 

7.71 

1.03 

fitl,all 

0.32 

7.65 

1.05 

Ht2,all 

1.04 

7.60 

1.00 

var,bath 

0.00 

7.76 

1.05 

var,all 

-0.30 

7.75 

1.02 


below the Fermi energy. Therefore only 22 (or 10 in the 
case of “bath”) parameters per orbital have to be varied 
to find the optimal unitary transformation to embed the 
effective model into the full Hilbert space. 

C. Results 

We will first compare free energy estimates as well as 
different local observables obtained from variational ED 
treatments to unrestricted Hartree-Fock as well as QMC 
calculations. Afterwards, we investigate the nature of the 
optimized effective bath and impurity states as obtained 
from the variational ED treatment. 


1. Free energy functional and local observables 

Tabic HIl shows the free energy functional $ (relative to 
unrestricted Hartree-Fock (UHF)), as obtained with dif¬ 
ferent starting points and different amounts of variational 
degrees of freedom in the variational ED approach. In 
the case of “bath”, the constant weight function (“fitO”, 
W n = 1) leads to the lowest values of the <F [p]. The mod¬ 
els derived using weight functions W n = l/w n (“fitl”) 
and W n = 1 /uj% (“fit2”) lead to free energy estimates 
which are about 1 to 3 eV higher in energy. The situ¬ 
ation for the case of “all” is similar. On this basis, we 
have chosen the starting guess obtained with the constant 
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weight function for the full optimization of the effective 
model parameters. The resulting parameters are shown 
in the last column of Tab. Q] and are close to the start¬ 
ing guess. The full optimization schemes (“var,bath” and 
“var,all”) find parameters which lower the functional $ 
considerably for “all” and slightly for “bath”. 

Regarding the impurity occupation ((rid), see Tab. HT1) . 
we see that the description by unrestricted Hartree-Fock 
is rather close to QMC, whereas restricted Hartree-Fock 
overestimates the occupation. All versions of exact di- 
agonalization lead to occupations close to the QMC re¬ 
sults, and many cases within the QMC error bars. The 
spin S (defined as (S' 2 ) = S(S + 1)), which is a two- 
particle observable, reveals the problems of the Hartree- 
Fock description. S is vastly overestimated by unre¬ 
stricted Hartree-Fock. The variational ED methods, es¬ 
pecially the “all” case for the constant weight function 
(“fitO”) and the full optimized ED model, lead to results 
close to QMC. 

To compare the results of the variational ED method 
with those ED methods based on fitting of the hybridiza¬ 
tion function on the imaginary axis, we should compare 
the “fitO,bath”,“fit 1,bath”, and “fit2,bath” cases to the 
corresponding “all” and “var,all” cases. We see that hav¬ 
ing more variational degrees of freedom leads to an im¬ 
proved description of the free energies, as it should be. 

In general, we learn that only in the case of optimizing 
both effective bath and effective impurity states (termed 
“all”) we reach lower values of the free energy functional 

than with unrestricted Hartree-Fock: The freedom to 
form mixtures of bath and impurity states in the effec¬ 
tive model is important to describe the free energy and 
local observables of the system adequately. As the vari¬ 
ational ED method provides more accurate (free) energy 
estimates than unrestricted Hartree-Fock, the approach 
introduced here could be a way to improve LDA+U total 
energy schemes. 


2. The effective basis states 


Now we analyze the unitary transformation relating 
the optimized basis states of the effective model and the 
original basis states. The transformation obtained for the 
models from the starting guesses with weight functions 
W n = 1, l/oo n , and 1/w 2 is shown in Fig. [3] b), c), 
and d), respectively. We observe a clear trend that the 
admixture of the original bath states into the effective 
bath state ( u , blue lines) is strongest in the vicinity of 
the effective bath site energy £ a \. For bath states close 
to the Fermi energy we get a sharp cut off for states on 
the opposite side of the Fermi energy. This is very similar 
to first order configuration interaction treatments of the 
AIM^. The original impurity admixture in the effective 
impurity (u d ) rises with the distance of the effective 
bath site from the Fermi energy. 




d) W n = 1/u>1 e *( eV ) 


c >2g £e g 
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Figure 3. (Color online) a) Hopping matrix elements between 
the impurity orbitals and the bath from the original AIM for 
Co impurities in Cu (solid t 2 g , dashed e g ). b)-d) Coefficients 
defining the optimal transformation from original bath states 
to effective bath states (blue/dark gray) and effective impu¬ 
rity states (red/light gray), c.f. Eqs. ([7J and (|8|). Optimized 
transformations for different effective models defined through 
fits of the hybridization with weight functions W n = 1 b), 
W n = 1 /uJn c), and W„ = 1/w 2 d) are shown. The ener¬ 
gies of the effective coupled bath sites i a are depicted as thin 
vertical lines. The numerical values of the transformation co¬ 
efficients defining the admixture of original impurity states to 
the effective bath and impurity states are given as insets. 


V. SPECTRAL FUNCTIONS 


The variational principle results in an effective model 
which represents thermodynamic ground state proper¬ 
ties in an optimal manner. This is a necessary but not 





































a sufficient condition to give a good approximation also 
for excitation spectra. In the following, we study the 
one particle spectral function for single orbital impu¬ 
rity benchmark systems from Sec. Mil with the bath 
states centered around = 0.3 eV and two different 
hybridization strengths, 14 = 0.9 eV and 14 = 0.3 eV, 
respectively. The impurity spectral function is obtained 
from the Lehmann representation of the impurity Green’s 
function 


G a (w) 


J_ I (/^Ma \ v ) I _ / -PE, 

Z ^ u> + E„- E.-i 0+ ^ 

fl V ^ 


(18) 


where in our calculations 0 + is replaced by a broadening 
of S = 0.1 eV and the inverse temperature is /3 = 3200, 
which is very close to the T = 0 calculations of expecta¬ 
tion values in Sec. IIIII 

We assess the quality of the spectra obtained from the 
variational ED method, from ED with the hybridization 
function fitted on the Matsubara axis (with two differ¬ 
ent weight functions W n = 1 and W n = 1 /u> n ) and from 
an unrestricted Hartree-Fock treatment by comparing to 
the exact spectrum of the original model. For the case 
of W n = l/oj n , we additionally compare the spectra for 
different amounts of variational freedom in choosing the 
basis states of the effective model, where we either op¬ 
timized the bath states only (termed “bath”, c.f. Sec. 
m or allowed impurity and bath states to mix (termed 
“all”). 

The dominant features of the original spectrum in the 
case of strong hybridization (14 = 0.9 eV, Fig. S]a)) are 
two major peaks at about —2.5 eV and 2 eV (stemming 
from bonding and anti-bonding combinations of impurity 
and bath orbitals), two satellite peaks far from the Fermi 
energy and additional smaller peaks around the Fermi 
energy. The spectral function from Hartree-Fock repro¬ 
duces the bonding/anti bonding peaks and those close to 
the Fermi energy very well, while the satellites are miss¬ 
ing. The variational ED method describes the positions 
of the main peaks well and also reproduces the satellite 
peaks, whereas the minor peaks around the Fermi energy 
are not present. The spectral function from a fit on the 
imaginary axis with a constant weight function (W n = 1) 
shows a similar picture but with major peaks and satel¬ 
lites shifted considerably towards the Fermi energy. The 
result for the weight function emphasizing small Matsub¬ 
ara frequencies ( W n = l/w„) leads to a good represen¬ 
tation of the peaks around —2.5 eV and 2 eV and some 
minor peaks around the Fermi energy but the satellites 
are missing completely. The resulting spectrum obtained 
by not mixing bath and impurity states (“bath”) shows 
only little resemblance to the original spectrum. I.e. in 
this case the mixing of bath and impurity basis states can 
not only improve total energies / thermodynamic poten¬ 
tials but also spectra quite significantly. 

In the case of weaker hybridization (14 = 0.3 eV, 
Fig. H]b)) the original impurity spectral function shows 


two Hubbard peaks at about —2.6 eV and 2.4 eV and 
in comparison to the former case more spectral weight 
and additional features close to the Fermi energy. Here, 
the Hartree-Fock description results in a spin-polarized 
ground state and describes the positions of the Hubbard 
peaks in the spectrum correctly. The enhanced spectral 
weight at the Fermi energy is however not reproduced. 
Exact diagonalization of the effective models obtained 
from the variational method and from the fits of the hy¬ 
bridization functions on imaginary axis give similar re¬ 
sults. In addition to the upper and lower Hubbard peaks 
the ED methods also reproduce enhanced spectral weight 
at the Fermi level. 

While the performance of the fit methods in reproduc¬ 
ing the spectra of the original model differs between the 
case with strong and weak hybridization particularly for 
the weight function W n = l/tu^, the variational method 
gives satisfactory results in both cases. The spectra from 
variational ED are in both cases at least as close to 
the spectra of the original model as the best spectrum 
obtained with any of the two bath fitting procedures 
W n = 1 or l/w n ) under investigation. 


a) 14 = 0.9 eV 



b) 14 = 0.3 eV w ( el/ ) 



Figure 4. (Color online) One particle impurity spectral func¬ 
tions for the single-orbital benchmark model introduced in 
Sec. uni with Eb = —0.3 eV and hybridization strengths 
14 = 0.9 eV (a) and 14 = 0.3 eV (b). The spectral func¬ 
tion from the exact solution of the original model is shown in 
bold cyan, from the Hartree-Fock calculation in red and from 
the variational ED method in black. Spectra from ED with 
fitted hybridization functions on the Matsubara axis are de¬ 
picted in dashed green (“fitO”, i.e. weight function W„ = 1), 
dotted blue (“fitl”, W„ = 1 /uj n optimizing all states) and 
dashed blue (“fitl”, W„ = 1/uj„ optimizing only bath states). 
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VI. CONCLUSION 

In conclusion, we present a variational exact diago- 
nalization method which provides self-consistently opti¬ 
mized parameters of discretized Anderson impurity mod¬ 
els. The method is based on an optimal partitioning of 
the system into a correlated part, where electronic inter¬ 
actions are explicitly taken into account, and an uncor¬ 
related rest and on finding optimal effective Hamiltoni¬ 
ans for both parts of the system. To this end, a varia¬ 
tion of a free energy functional w.r.t. one-particle basis 
states spanning the correlated subspace and the matrix 
elements of the effective Hamiltonians is performed. 

A benchmark of the variational ED method against an 
exact solution of a one orbital Anderson model demon¬ 
strates its excellent performance in reproducing ground 
state observables of the impurity and bath and addi¬ 
tionally a sound performance in reproducing the im¬ 
purity spectral function. A comparison with Hartree- 
Fock and established bath discretization schemes for ED 
shows that the variational approach introduced here even 
works for difficult cases, i.e. when the bath is symmet¬ 
ric around the Fermi energy. Furthermore, applicability 
of the variational ED method to realistic multi-orbital 
cases is demonstrated with the example of Co impurities 
in bulk Cu. Also here, the variational method leads to an 
accurate description of local one and two particle observ¬ 
ables like the impurity occupation and the spin. Energet¬ 
ically the method outperforms unrestricted Hartree-Fock, 
which suggests that the variational ED approach could 
be useful to improve total-energy approaches to corre¬ 
lated systems beyond LDA+U. Finally, the method intro¬ 
duced, here, can be used to embed established bath dis¬ 
cretization schemes such as the fit of hybridization func¬ 
tions on Matsubara frequencies into a variational frame¬ 
work and to reach an unbiased decision of e.g. which 
weight function to choose in the fits of the hybridiza¬ 
tion functions. In the example studied, here, the con¬ 
stant weight function leads to best results in terms of the 
free energy, whereas W n = 1/w 2 leads to results quali¬ 
tatively similar to a first order configuration interaction 
expansion^. For the description of systems closer to the 
atomic limit, like Co on Cu or 4f systems we expect Cl 
expansions to converge faster and weight functions like 
W n = 1 /w n or W n = 1/w 2 might be a better choice. The 
presented variational approach will allow for an unbiased 
decision in any case. 

Additionally, the variational method is universal and 
different choices of the effective correlated space are pos¬ 
sible. To gain yet higher accuracy, more bath levels could 


be included. On the other hand, to make contact with 
ligand field theory or crystal field theory descriptions of 
magnetic impurity systems and nanostructures one could 
consider correlated subspaces with very few or without 
any effective bath orbitals at all. 
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Appendix A: Fit of hybridization functions on the 
Matsubara axis 

The method of fitting hybridization functions is shortly 
introduced for the sake of completeness. The basic idea 
is to minimize a cost function for the inverse impurity 
Green’s function (or equivalently the hybridization func¬ 
tion) of the discretized model and that of the original 
model, both defined on the imaginary frequency axisi£. 
In the case of one effective bath site, the discrete impurity 
Green’s function is defined as 

( u 2 V 1 

g 0 (iuJn) = - e d - n - r-— (Al) 

y iuj n - £i J 

and the Green’s function of the original model as 

) -i 

■ ( A2 ) 

The cost function then reads 

-t Umax 

x 2 =-—T 5Z Wn l^o - go\iu n )\~, (A3) 

rimax i -L 

n —0 

where W n is a weight function. Popular choices for 
the weight function are W n = 1 , W n = l/w n and 
W n = 1/w 2 . Different weight functions put different em¬ 
phasis of low/higher Matsubara frequencies^. Through¬ 
out this work, we have chosen /3 = 40 and n ma x = 1000. 
This method only provides the effective parameters i\ 
and V. However, in order to calculate the functional 
<f>[jO^] an optimal unitary transformation in above sense 
is calculated and the hj^ k , are found by a projection of 
a Hartree-Fock solution onto the basis states of 1Z. We 
assume that the effective energy of the impurity site is 
the same as in the original model {id = £d). 
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